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ABSTRACT. The next generation ground-based telescopes rely heavily on adap- 
tive optics for overcoming the limitation of atmospheric turbulence. In the future 
adaptive optics modalities, like multi-conjugate adaptive optics (MCAO), at- 
mospheric tomography is the major mathematical and computational challenge. 
In this severely ill-posed problem a fast and stable reconstruction algorithm is 
needed that can take into account many real-life phenomena of telescope imag- 
ing. We introduce a novel reconstruction method for the atmospheric tomog- 
raphy problem and demonstrate its performance and flexibility in the context 
of MCAO. Our method is based on using locality properties of compactly sup- 
ported wavelets, both in the spatial and frequency domain. The reconstruction in 
the atmospheric tomography problem is obtained by solving the Bayesian MAP 
estimator with a conjugate gradient based algorithm. An accelerated algorithm 
with preconditioning is also introduced. Numerical performance is demonstrated 
on the official end-to-end simulation tool OCTOPUS of European Southern Ob- 
servatory. 



1. INTRODUCTION 

The resolution of an optical imaging system can be limited by several factors. 
Whereas imperfections of the optical setting can be improved, diffraction always 
defines a fundamental limit in this regard. In ground-based telescope imaging 
with a small mirror, diffraction is often the dominating effect. However, the in- 
fluence of the atmospheric turbulence scales much faster with the increase of the 
mirror diameter. Already in the next generation of telescopes, called the extremely 
large telescopes, atmospheric turbulence is the major limiting factor for the angu- 
lar resolution far beyond the diffraction-limit. It is a great challenge for science 
and technology to find ways to achieve diffraction-limited imaging for the future 
ground-based telescopes. 

During the last decades adaptive optics (AO) technology has developed into a 
powerful remedy for this problem. Adaptive optics refers to real-time compensa- 
tion for the distortions in the wavefronts of incoming light due to the atmospheric 
turbulence. Although the technology involves a number of engineering challenges, 
the benefits have proven to be fundamental. Adaptive optics correction has been 
implemented in many major telescope projects, e.g., the Very Large Telescope and 
the Gran Telescopio Canarias. Moreover, AO is planned as an essential part of all 
extremely large telescopes. The work for this paper was largely carried out within a 
project established towards developing mathematical algorithms for the European 
Extremely Large Telescope (E-ELT), the future telescope of the European Southern 
Observatory (ESO). 

The mathematical challenges and future prospects for inverse problems field 
were comprehensively reviewed by Ellerbroek and Vogel in lfl4l . With the arrival 

1 



2 



TAPIO HELIN AND MYKHAYLO YUDYTSKIY 



of next generation implementations of adaptive optics, the severely ill-posed atmo- 
spheric tomography problem becomes the crux of mathematical AO research. In 
this paper we discuss atmospheric tomography in the context of multi-conjugate 
adaptive optics (MCAO). 

In the classical AO, a single guide star, i.e., a single light source is observed. 
The aberrations in the incoming light are measured which enables the AO system 
to correct the cumulative effect of the turbulence towards directions close to the 
guide star. The correction is performed with a deformable mirror (DM). Multi- 
conjugate adaptive optics extends this idea by using several guide stars and mul- 
tiple deformable mirrors. First, the data obtained from the guide stars is used to 
solve a reconstruction of the turbulence profile in the atmosphere (see Fig. [I]). This 
problem is called atmospheric tomography. Second, having multiple deformable 
mirrors and the three dimensional reconstruction of the turbulence enable the as- 
tronomer to correct for much larger field of view than in the classical implementa- 
tions. 

The physics of turbulence is an extensive field of study and much is understood 
about how the turbulence in the atmosphere is formed. Statistical models for tur- 
bulence are frequently utilized in the adaptive optics literature by postulating the 
tomography step as a Bayesian inference problem. In addition, this makes it possi- 
ble to take into account the statistical nature of the measurement noise. The max- 
imum a posteriori (MAP) estimate is the standard point estimate used to describe 
the resulting Gaussian posteriori distribution. 

In the MCAO related literature, both iterative and non-iterative solution meth- 
ods have been proposed for the MAP estimator. With the launch of the next gener- 
ation of extremely large telescopes the dimension of the problem increases rapidly. 
Even with heavy parallelization, the non-iterative methods have a high compu- 
tational cost and the research in recent years has been inclined towards iterative 
methods. 

The iterative solvers of the MAP estimate are typically based on conjugated gra- 
dient (CG) methods, see [Q2] and references therein. Effort has been put into devel- 
oping an efficient preconditioner for the problem. The multigrid preconditioners 
are investigated in Il9ll20l[l8l . whereas in B4ll45l . Fourier domain precondition- 
ers have been proposed. 

Especially for iterative methods it is of value to be able to represent the opera- 
tors in a sparse form. In this regard the Fourier basis is very useful and Fourier- 
transform based reconstructors have been proposed in B3ll42l[T6ll . In other typical 
bases, the forward operator and the inverse covariance of the noise are fast to ap- 
ply. However, the inverse covariance of the prior is often a full matrix and thus a 
sparse approximation is required. In the approach introduced by Ellerbroek ifTTTl 
the turbulence power law is modified in order to achieve a sparse approximation by 
biharmonics. Later, a very promising CG based method called the Fractal Iterative 
Method (FrIM) has been developed by Tallon and others in [[39l |4H [381 El - There, 
the inverse covariance is approximated by a factorization that can be applied in 
0(n) operations. 

We also point out that outside the Bayesian framework the atmospheric tomog- 
raphy problem has been approached by an algorithm based on the Kaczmarz itera- 
tion. The method was introduced by Ramlau and Rosensteiner in ll32l . The authors 
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obtain a very efficient matrix-free solver by performing the wavefront reconstruc- 
tion and the tomography step in two separate problems. The reconstruction of 
the atmosphere from incoming wavefronts by the Kaczmarz iteration delivers very 
promising results in good imaging conditions. The incoming wavefronts are recon- 
structed from the measurement with an algorithm called the CuReD ll35l . Several 
solvers for the wavefront reconstruction exist (see, e.g., Il2l l30l '). 

In this paper we suggest a method utilizing compactly supported orthonormal 
wavelets to represent the atmosphere. Our method is a CG based iterative method 
that solves the MAP estimate for the atmospheric tomography problem. We demon- 
strate a successful performance in low flux imaging conditions, i.e., when only 
a low number of photons can be measured, and with respect to some practical 
phenomena that are well-known to limit the reconstruction quality. We have im- 
plemented the method with the Daubechies wavelets [ 8 , 6 ] in order to have good 
reconstruction of local details. In addition, the Daubechies wavelets have a use- 
ful localization in the frequency domain. This enables us to introduce our key 
contribution by approximating the inverse covariance with completely diagonal 
representation. It is shown that such a representation produces an equivalent reg- 
ularization term in the MAP estimation problem as indicated by the theory. We 
point out that this approximation is flexible with respect to choosing a different 
model for the turbulence power law. In terms of temporal control we rely on an 
established method called the pseudo-open loop control (POLC) which has been 
demonstrated to be very robust [29]. 

In the numerical tests we introduce two variants of the algorithm. In the first 
setting, we reconstruct more layers of the atmosphere than deformable mirrors us- 
ing the CG algorithm and optimize the DM shapes accordingly. In our opinion this 
demonstrates well the best qualitative performance of the wavelet based method. 
Moreover, we investigate the stability of the regularization procedure in this setting. 
The second variant of the method is an accelerated algorithm developed towards 
achieving the real-time requirements. Here, layers are reconstructed at the alti- 
tudes of the deformable mirrors and DM shapes are chosen as the reconstructed 
layers. In the accelerated method we utilize a modified Jacobi preconditioner, for 
which we demonstrate fast convergence. Numerical simulations are carried out 
on the OCTOPUS, the official end-to-end simulation tool of ESO. We illustrate 
the performance of our method in the low flux imaging conditions and compare 
these results against the matrix-vector multiply (MVM) algorithm, which is the 
benchmark reconstructor of ESO. 

Wavelet methods in adaptive optics have been previously studied in ll22l |2D . 
however not in the context of MCAO or in the Bayesian scheme. In the field of 
inverse problems wavelets are applied widely (e.g., ll37l l24lD . For an extensive 
introduction to wavelet basis we refer to |9|. 

Notice that atmospheric tomography is a severely ill-posed problem and is very 
closely connected to limited angle tomography iflOl . Thus, from the general per- 
spective of inverse problems, the theoretical limitations of MCAO are interesting 
and have been considered in ||43ll42l . Inverse problems related to waves travelling 
in random media have been considered, e.g., in the works of Papanicolaou, Bal and 
Borcea S LULLS. 
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This paper is structured as follows. In Section [2] we discuss the mathematical 
model for the propagation of light through the atmosphere and how the measure- 
ments are obtained. We close the section by explaining the Bayesian paradigm and 
how the MAP estimate for the atmospheric tomography is solved. In Section [3] 
we introduce the diagonal approximation for the inverse covariance of the turbu- 
lence statistics. Section [4] features parts of our method that are essential to MCAO 
solver, but which have been studied before. Here, we discuss the fitting step and 
the control algorithm. The concepts of spot elongation and tip-tilt uncertainty are 
also introduced. These practical phenomena have an essential impact on the noise 
model. Finally, in Section [5] we demonstrate the numerical performance of our 
method. 



2. Atmospheric tomography 

2.1. Problem setting. The wind in the atmosphere causes an irregular mixing of 
warm and cold air. This effect is called the atmospheric turbulence. The fluctu- 
ations of the temperature are essentially proportional to the refractive index fluc- 
tuations ll33ll and hence the turbulence affects the propagation of light. With the 
geometric optics approximation and under appropriate assumptions on the atmo- 
sphere, the phase of light at the aperture is distorted according to 

(1) #r,0)« I" p(r + 6 ■ H)dH, 

Jo 

to a good approximation towards directions 6 = (61,62, 1) close to the zenith. 
Above, p describes the fluctuations of the refractive index, r = (x, y, 0) is the 
location at the aperture, and H is the height of the atmosphere. The approximation 
([I]) is derived in Il40ll36l . The challenge in atmospheric tomography is to obtain a 
good estimate of p based on indirect measurements of (f>(-, 6) towards directions 
of the guide stars. 

The strength of the turbulence at a given altitude varies heavily. However, at 
the typical telescope sites most of the turbulence is concentrated on certain alti- 
tudes. This observation has given rise to a layered atmospheric model, where the 
refractive index is approximated on a finite number of two-dimensional layers at 
fixed altitudes. In the simplest example, in ground layer adaptive optics only one 
layer is considered since the majority of the turbulence strength is located close to 
the aperture of the telescope (the ground layer). Due to the availability of several 
deformable mirrors, the implementation of MCAO benefits from a more accurate 
description of the atmosphere. 

Let us consider how the equation ([T]) reduces for a layered atmosphere model. 
We denote each modelled layer, located at altitude hi by cj)£, £ = 1, . . . , L, and 
by <f> = (</>i> • • • 5 <^l) a vector representing the atmosphere. Assuming geometric 
propagation, the light arriving from infinity produces incoming wavefronts accord- 
ing to 

L 

(2) 0(r, e) = P NGS = Y, Mr + Ohe), 

where r denotes a point inside the aperture and the vector 6 describes the direction 
of the guide star (see Fig. Mb. Here, the projection Pq GS maps the atmosphere to 
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the incoming wavefront from direction 6. More details on the wave propagation 
through the layered atmosphere model can be found in [34]. 

However, in practice, there are not enough bright natural guide stars to cover 
most areas of interest. To overcome this problem, astronomers have developed 
a technology utilizing lasers, which can generate artificial stars at finite altitude. 
A laser guide star (LGS) is produced by shooting a powerful laser into the atmo- 
sphere where it scatters strongly at certain higher altitudes (33). Due to the finite 
altitude, the light arriving to the telescope passes through a cone-like volume in 
the atmosphere (see Fig. [2]>. The corresponding distortion then satisfies 

^(r, e) = p LGS = fa ((i - ^) r + ^) » 

where H denotes the altitude where the laser scatters. Again, Pg GS stands for the 
projection of the atmosphere to the incoming wavefront with respect to the cone 
geometry. Whereas an LGS provides a very bright source at directions where no 
NGS is available, it also introduces some practical limitations. These phenomena, 
called tip-tilt effect and spot elongation, are described in Section|4] 

The incoming wavefront can be measured indirectly. A common measurement 
device in adaptive optics is the Shack-Hartmann wavefront sensor, which was de- 
scribed in detail in [14]. Essentially, a Shack-Hartmann sensor measures a quantity 
proportional to the average gradient of the wavefront on a rectangular grid formed 
by small lenslets, i.e., 



(3) 8ij = C V0(r, 0)dv 

where sa = (sf,-,sf-) G M? is the measurement and C is a constant. Above, 
Qij denotes the lenslet also referred to as the subaperture. Other wavefront sensor 
modalities exist (e.g., curvature sensor 11331 ) but are not considered here. 

Let us next describe an equation connecting the measurements with the atmo- 
sphere for a full MCAO system. The MCAO system we consider will utilize both 
laser and natural guide stars. Moreover, each guide star here has an individual di- 
rection which we often use to index both the guide star and their corresponding 
wavefront sensor (WFS). Next, we use T to denote the measurement operator 

(4) S0 = rvH-,0) 

where is a vector containing all grid values Sij from formula ([3]). The Shack- 
Hartmann sensors modelled in equation Q can have different resolution and hence 
Tg is direction-dependent. 

In what follows we consider a system that observes M guide stars. Their di- 
rections are denoted by 6 m , where 1 < m < M. We assume that out of this 
number, the first Mlgs are laser guide stars. For the rest of the paper we simplify 
the direction-dependent notations by replacing 6 m by m whenever no confusion 
appears. Now we can write the subproblems for different guide star directions by 

(5) s m = r m pLGS0 and Sm , = r m ,pNGs^ 

for 1 < m < Mlgs an d Mlgs < m' < M. The full system is then described by 

(6) s = (s m )^ = i = A0, 
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Figure 1 . In atmospheric tomography the turbulence layers are 
reconstructed from the wavefront sensor data. 
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Figure 2. The cone effect. Laser guide star light source is fixed 
above the turbulence layers at some finite altitude. Light traveling 
from the LGS to the telescope pupil passes through smaller areas 
at higher turbulence layers. 



where A is the concatenation of operators T m P^ s and T m /P ? ^P s . Estimating 4> 
from a given s is called the atmospheric tomography problem. 

2.2. Bayesian inference. The Bayesian inference is a standard approach for solv- 
ing problem ([6]). This appears natural since information is available about the 
statistical behavior of the unknown wavefronts and the measurement noise. The 
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Bayesian paradigm considers problem ([6]) in a random setting S = A<J> + £, where 
S, <& and £ denote random variables describing the measurement, the incoming 
turbulence wavefront and the additive noise, respectively. Given a sample of S, 
i.e., the measurement, the task is to deduce information about the unknown <J>. 
Both $ and £ are typically modelled as Gaussian random variables. We discuss 
the distribution of the wavefronts and their physical interpretation below in Section 
[3] in more detail. The noise in the Shack-Hartmann measurements is produced by 
several components ll33l . However, for the LGS measurements, the effect of spot 
elongation has a major influence on the noise distribution. We return to the noise 
covariance and the spot elongation in Section [4] 

Below, we denote the covariance operators of <J> and £ by C$ and Cg, respec- 
tively. Furthermore, it is assumed that the separate layers are zero centered and 
uncorrected. This implies that C$ has a block-diagonal structure 

C$ =diag(Ci,...,Cx), 

where Q denotes the covariance of the layer I. In the setup given above, the 
maximum a posteriori estimate can be obtained by solving 

(7) = argmin (||C^ 1/2 ^||| + ||CJ 1/2 (s - A0)|||) . 
For a general introduction to the Bayesian inverse problems, see IT231 . 

3. Turbulence models and the Bernstein-Jackson equivalence 

Let us consider for the moment the theory of Gaussian random variables in real 
separable Hilbert spaces. Let <j> be a measurable map from the probability space 
f2 to a Hilbert space H. Then is Gaussian if and only if for all p\ , . . . , p m G H 
the mapping 3 uj \- > ({4>, pj))™ =1 is a Gaussian random variable in M m . The 
distribution of <fi is determined by the expectation E(p and the covariance operator 
: H — » H defined by 

WuCfih) = E ((0 — Ecj>, ^i)<0 - Ecf>, ^)) . 

It is well-known that for any ip € H and a linear positive self-adjoint trace class 
operator C in H with M{C) = {0} there exists a Gaussian random variable cj) in 
H with mean tp and covariance C OTTl . 

Below, we are concerned with zero-centered random variables <fi that have real- 
izations in some Sobolev space H S (M?) with s > and a covariance operator of 
the form 

(8) C^ = F*MF. 

Above, F is the Fourier transform on M 2 and M is a multiplication operator 
Mf{n) = m{n)f{K) where m is a positive bounded function. With an appropriate 
decay of m at infinity, the operator is trace class and hence <fi is well-defined. 
Further, the Schwartz kernel k^x, y) of the operator satisfies 

k^(x, y) = E (<p(x) - E<f>(x)) {<t>{y) - Et/>(y)) 

in the sense of generalized functions. We call the covariance function of the 
random process <fi. 
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In the literature of adaptive optics, a turbulence layer <fi is assumed to be isotr opic 
and stationary, i.e., depends only on the distance of x and y. In particular, is 
completely characterized by 

(9) k<f,(z) = k ( f ) (x,x + z). 

where x, z £ M 2 . Now, it can be shown that if the statistics of 4> satisfy equations 
® and ©, then 

m(/c) = (Fk^)^) 

in the sense of tempered distributions. The multiplier function m is often referred 
to as the power spectrum. The classical model based on the Kolmogorov-Obukhov 
law of turbulence states that the power spectrum m follows a power law 

(10) m(«) = C\k\- 11/3 

inside the so-called inertial range Ki n < \n\ < K out with some constant C. It is 



not straightforward to expand the power law ( pL0] > outside the inertial range due to 
the strong singularity at zero. We make a common choice of von Karman model 
(33l[T4l to modify (TO]) by assuming 



1 \ - n /6 

1 ■ \2 



m(n) = c p (h) I —2 — h \k\ 

\ K out 

where K ou t is the outer scale of the turbulence and c p (h) describes the measure of 
the optical turbulence strength depending on the altitude. This choice for the power 
law coincides asymptotically with ( fTO] ) in the high-frequency regime, however, the 
singularity at zero is removed. In conclusion, we notice that an equivalence 



(ii) \\c^ /A f\\h = WiKut + N 2 )*J7ll£» * KJ WfWb + M-A)*f\\h 

holds in the Cameron-Martin space of <f>, i.e., for any / £ ^ 11//6 (M 2 ). Here and in 
what follows, we write p ~ q if the two pseudo-norms p and q are equivalent. 

Assume that the wavelets studied here are r-regular, i.e., have r vanishing mo- 
ments and are r times continuously differentiable. For sufficiently large r the last 



term in ( 1 1 1 is equivalent with the expression 



(12) ||(-A)I2/|| 2 - 2 ^^2 2 -^|(/,Va)| 2 , 

A=l 

where j is the wavelet scale of the wavelet ip\ with global index A. The equivalence 
above is known as the Bernstein- Jackson inequalities [27]. 

In the discretized problem, the function / in ( [TT| ) is represented by finite number 
of wavelet scales. In that case, an equivalent representation for the regularizing 



term in ( 1 1 1 can be produced by a diagonal matrix Di : M. ne — > M ra< that satisfies 



1 



ii 



ii 



Above, ri£ denotes the total number of wavelets for layer I. Moreover, we denote 
D = diag(Z?i, . . . , Dl). Finally, by approximating the prior covariance of the 
discretized problem by the ideal model we get 

L L 

(13) ||C $ 1/2 0II( L2) l = Yl \\C; 1/2 Uh - E(^ c ^) 2 = ( Dc ' c ) 2 
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for 4> = {<t>e)e=i an d the wavelet decomposition = Y^x=i c ^,aV^,A with respect 
to the wavelet basis {ip£\ '■ A = 1, . . . ,nf\ of layer t. Above, we denote by 
c £ = (q,a)aLi the vector of wavelet coefficients associated to layer £, and by c the 
concatenation of vectors eg. 

— 11/3 

We point out that in practice the term K out becomes negligible. Furthermore, 



the approximation error introduced in < [T3[ > is beyond the scope of this paper. In 
numerical simulations presented in Section [5] we study how different weighting of 
the regularizing term (Dc, 0)2 affects the reconstructions obtained by the method. 

4. Other features of MCAO 

4. 1. Fitting step. In an MCAO system the correction for the wavefront distortions 
is produced by several deformable mirrors that are conjugated to different altitudes. 
Hence, a successful mathematical algorithm for MCAO must also deduce optimal 
mirror shapes for the DMs based on the reconstruction from equation ([6]). This 
subproblem is called the fitting step. Given a sufficient reconstruction of the at- 
mosphere, the fitting step is a well-posed least squares minimization problem and 
thus classical solution methods provide an efficient reconstruction strategy. We 
point out that in the ideal fitting one aims to minimize a functional 

(14) argmin a E(/ / (H g a - Pg GS (f>) 2 dxdGj , 



f Jn 



where a is the correction profile, Hg is the collection towards direction 6 and Pg GS 
is defined in equation ([2]). Moreover, Q is the aperture domain and 6 belongs to 
the field of view F. The problem is typically discretized by choosing a finite set of 



directions over which the difference in equation ( |T4| ) is averaged. We follow this 
tradition by formulating the fitting step as the minimum norm solution to 

(15) argmin a ||Ha — P0|| 2 

where H and P are concatenations of operators Hj and Pj 1 " 35 , respectively, to- 
wards a finite set of directions Oj, j = 1, . . . , N sampled from the field of view. 
For more detailed discussion on the fitting step see lfl4l . 

4.2. Closed loop control. Although in next generation telescopes the DMs are 
adjusted within milliseconds, the delay between the measurement and the applied 
DM correction induces an error that needs to be considered. Consequently, a robust 
temporal control is a fundamental part of the system. 

In an MCAO system, the wavefront sensors are located behind the deformable 
mirrors in the optical path of light. This is contrary to our assumption on the prior 
model discussed in the previous section, as the WFS measures the residuals of the 
incoming wavefronts, instead of the incoming wavefronts themselves. This mode 
of operation is called closed loop. 

In order to model the physics of turbulence in the prior covariance, we follow 
here a method called the pseudo-open loop control (POLC). The straightforward 
idea of POLC is to approximate open loop measurements by combing the mirror 
shapes with the residual data. The POLC was introduced to AO in lfT3l and further 
studied in ifTTl [29|| . It has proven to be stable and robust against large levels of 
system errors [29]. 

We rely on a modified POLC, where an integrator is used in the control scheme 
(see e.g., [29]). We assume that our system has a two time-step delay. Let t G 
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N, t > 2 denote time-steps. Further, residual Shack-Hartmann data s are measured 
over the time period [t — 2, t — 1) for the mirror corrections a t _2- We assume that 
the reconstruction step (computing DM shapes from measurements) consumes one 
time period, [t — l,t) and the computed mirror shapes a t are applied to the mirrors 
at time step t. 

Algorithm 4.1. Pseudo-open loop control. 

(1) s 01 = s + Aa t _ 2 

(2) Aa = Rs o1 - a t _ 2 

(3) At = a t _i + gAa 

Above, A maps the mirror shapes a to the correction in the measurement space, 
similar to Moreover, the reconstruction operator R maps WFS measurements 
to DM shapes. The scalar parameter < g < 1 denotes the gain, which controls 
the mirror update. 

We point out that the POLC is not optimal in the sense of cumulative residual 
variance. More sophisticated Kalman filter based methods can achieve this [28]. 
The drawback however in such methods is the additional computational load that 
is a limiting factor, especially for MCAO. Numerically efficient solutions towards 
this end are an interesting topic for future research. 

4.3. Spot elongation and tip-tilt uncertainty. Measurements observed from an 
artificially generated laser guide star suffer from two deficiencies: spot elongation 
and tip-tilt uncertainty. These effects were discussed in detail in |[T4l and we only 
briefly state how we correct for these errors in our algorithm. A successful mathe- 
matical model must take these effects into account, as the performance of the AO 
system would degrade otherwise. 

The spot elongation effect occurs due to the physics of scattering at the sodium 
layer. In practice, a laser guide star is not an ideal point source but rather an 
extended three dimensional source. This translates to an elongation of the observed 
spot on the measurement device, which can be described by a correlation of x- and 
y-measurements in each individual subaperture of the Shack-Hartmann WFS. Our 
method handles spot elongation following the approach taken in . 

Let us give a brief overview of the noise covariance matrix Cf . For an MCAO 
system, which relies on a combination of laser guide star and natural guide star 
wavefront sensors, the full noise covariance matrix is given as a block-diagonal 
matrix with respect to the sensors, 



Here we associate a noise covariance matrix C m for each direction of the guide 
stars 6 m , m = 1, . . . , M. Recall that the first A/lgs directions are associated with 
sensors observing laser guide stars. The remaining natural guide star sensors are 
not affected by spot elongation. For those directions we assume that the noise is 
identically and independently distributed in all subapertures with variance a 2 . 

Also, the spot-elongated measurements are uncorrected between different sub- 
apertures. However, for any subaperture the noise in the x- and y-measurements 
is correlated and hence the covariance matrix is block diagonal, composed of 2x2 
blocks 
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where S is the total number of subapertures of the WFS in direction 6 m . Each 
block C m: i, i = 1, . . . , S, can be expressed as 

(16) C m:i = a 2 (i + j^P ■ P T ) 

where / is the full width at half maximum (FWHM) of the non-elongated spots, (3 
is the elongation vector and a 2 is like above. Moreover, the parameter < r < 1 
is used to tune the relative increase of the noise with respect to the elongation. 
In consequence, the block-diagonal structure of the covariance matrix implies that 
applying or its inverse is computationally very cheap. For details on deriving 
C m .i we refer to Q and the references therein. 

The tip-tilt uncertainty is closely related to the uncertainty of the location where 
the LGS scatters in the atmosphere. Such an error has a large impact on the wave- 
front sensor measurements. However, it can be shown that the major part of the 
error is contained in the average x- and y-derivatives over the whole sensor, i.e., 
the tip and the tilt of the incoming wavefront. 

There are several tip-tilt correction methods that can be applied, such as a split 
tomography approach fl8l . a coupled-equation approach lfl4l or a noise-weighted 
approach [44J . We use a more straightforward approach, in which we remove the 
incorrect tip-tilt component in the LGS measurements by modifying equation ([5]) 
to 

(17) (I-T) Sm = (I-T)T m P^ GS (j> 

for m = 1, . . . , Mlgs> where T is an orthonormal projection into the tip and tilt 
components. Other way of stating ( fT7| ) is to say that we use 

(18) C- l = (I-T)C~ l {I-T) 

as the covariance matrix instead of C7" 1 . Hence this approach neglects more in- 
formation as, e.g., the noise-weighted approach and relies more on the NGS mea- 
surements. The successful performance of this method is supported by numerical 
tests. 

5. Numerical implementation 

5.1. Simulation environment and algorithm. For the simulations, we use the 
proposed multi-conjugate adaptive optics configuration for the European Extremely 
Large Telescope. The telescope gathers light through a circular pupil with diame- 
ter of 42 m, of which roughly 28 percent are obstructed. There are six laser guide 
stars positioned in a circle with a diameter of 2 arcmins. To each laser guide star, 
a Shack-Hartmann wavefront sensor with 84 x 84 subapertures is assigned. More- 
over, there are three natural guide stars positioned in a circle with a diameter of 
8/3 arcmins. The sensors assigned to the natural guide stars are low resolution 
Shack-Hartmann sensors (one with 2x2 and two with lxl subapertures) for tip- 
tilt correction. Further, the E-ELT uses a configuration of three deformable mirrors, 
located at altitudes km, 4 km and 12.7 km. The mirrors are modeled by piecewise 
bilinear functions, with a total number of 9296 degrees of freedom. 

We demonstrate the performance of our method on OCTOPUS ll25l . the official 
end-to-end simulation tool of the European Southern Observatory. The software 
generates nine frozen layers of the atmosphere located at altitudes between 47 m 
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Figure 3 . Star asterism of the six laser guide stars and three nat- 
ural guide stars, as well as the 5x5 quality evaluation grid over the 
field of view (in arcsec). 



and 1 8000 m. The evolution of the atmosphere is simulated by shifting these layers 
according to their wind directions and speed. 

In the test cases below we simulate one second of evolving atmosphere. The 
Shack-Hartmann measurements are read out 500 times per second. A two-step 



delay is observed, as described in Section 4.2 The measurements suffer from 
photon noise and readout noise. The quality of the reconstruction is evaluated in 
25 directions arranged in a 5 x 5 grid over the field of view, which is a square of 2 
arcmins. As a criteria the long exposure (LE) Strehl ratio 11331 is used in K band 
(for a wavelength of 2200 nm). The Strehl ratio is a commonly used measure of 
quality in the astronomical community. Towards directions close to the zenith it 
can be estimated by the Marechal approximation [33] according to 



s(0) « e 



-(2^(0)/A) 2 



where s is the Strehl ratio, r(0) is the root mean square error in the correction 
of the incoming wavefront from direction 6, and A is the wavelength. The long 
exposure Strehl relates to the average Strehl ratio over the observed timespan. The 
star asterism, as well as the 25 evaluation directions are depicted in Fig. [3] 

We base our algorithm on the POLC method described in Section |4~2{ where the 
operator R combines the solution operator to the atmospheric tomography prob- 
lem ([7]) and to the fitting step equation (15 1. In the special case of reconstructing 



turbulence layers directly at DM altitudes, we omit the fitting step equation and 
determine the mirror shapes by the shape of the reconstructed turbulence layers. 

The atmospheric tomography problem ([7]) discretized in the wavelet basis is 
equivalent to solving the linear system of equations 



(19) (A T C £ 1 A + aD)c = A T C £ r 1 



s. 



where A is the discretization of (roj). The role of the scalar parameter a is further 



discussed in Section pT2] We solve equation ( [19] ) using either the conjugate gradi- 
ent (CG) method or the preconditioned conjugate gradient (PCG) method with a 
modified Jacobi preconditioner, discussed below. 
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In the numerical simulations we utilize the Daubechies wavelet basis. They are 
a well-known orthogonal wavelet family with compact support (H and have a 
useful time-frequency localization. For a sufficiently large n, the Daubechies n 
wavelets are 2-regular and fulfill the equivalence ( fT2] ). In order to enhance the 
spatial localization we have chosen to use n = 3. It is well-known that the 
Daubechies 3 wavelets belong to the Holder space C 1+5 (R 2 ) with 6 « 0.0878 @. 
Even though they are not 2-regular, they form a Riesz basis in H 2 (M. 2 ) Q. Based 
on our numerical tests we believe that this is sufficient in practice. 

By formulating the problem in the wavelet domain we gain a significant im- 
provement in terms of convergence of the CG algorithm. This follows from the 
underlying operators possessing a favorable spectral structure. Further, the conver- 
gence is accelerated by using a modified Jacobi preconditioner, which we discuss 



in the following. The operator appearing on the left-hand side of equation ( 19 1 is 
given by 

M LGS M 

El TLGSnT^-1 TLGS i / TNGSnTPy— 1 TNGS i „t-v 

\ A m ) °m An "+" /—i ^ m ' ' m' rn' "+" aJJ ' 

m=l m'=M L GS+l 



where C^ 1 is defined by equation ( p~8] >, 



A^ S = r m pLGS^-l and A NGS = Tm p^S w -l 

Above, W^ 1 is the inverse wavelet transform mapping wavelets to functions and T m 
is the discretization of the Shack-Hartmann operator according to the Fried geom- 
etry (see, e.g., [ 33 ]). Following the discussion of Ellerbroek and Vogel in lfl4ll . we 
choose our preconditioner based on only the LGS components, as the low-rank 
perturbations, corresponding to only a finite number of eigenvalues, do not affect 
the asymptotic convergence rate of the conjugate gradient algorithm. Thus, our 
modified Jacobi preconditioner is 

/Mlgs \ 
J = diag ($£ S ) T C- l A h ™ + «D. 

\m=l / 

Finally, we reduce the number of conjugate gradient iterations by choosing the 
initial guess as the reconstruction in the previous time-step. This widely used 
technique for iterative methods in adaptive optics lfl4ll is known as warm restart. 

5.2. Stability of the regularization. The diagonal regularization operator in our 
method was obtained by using the Bernstein-Jackson equivalence. Clearly, the 
argument applied here does not state explicitly which value for a in ( fT9] ) is the 
optimal choice. Also, in this context a can be seen as a regularization parameter. 
Increasing its value can be considered as stabilization against modeling errors. We 
point out that there are several components of the problem that affect the stability, 
e.g., the temporal control (gain) and modeling of spot elongation. However, too 
large value will reduce the quality of the reconstructions. In the following, we 
demonstrate the performance of the method when a varies. 

We study a realistic noise-contaminated situation where the LGSs illuminate 
100 photons per subaperture and time-step. Furthermore, the spot elongation and 
tip-tilt uncertainty for the LGS are simulated. The NGS tip-tilt sensors observe 
500 photons per subaperture and frame. The readout noise is set to 3 and 5 electrons 
per pixel for the LGS and NGS sensors, respectively. 
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Figure 4. Long exposure Strehl vs. varying a in equation ( 19 ) 



Our algorithm is set as follows. We reconstruct nine layers at the altitudes of 
the nine simulated turbulence layers using the CG algorithm with 10 iterations. 



Further, we solve the fitting step equation (15 1 using the CG algorithm with 4 
iterations for optimization directions given by the 5x5 evaluation grid. We choose 
a gain of g = 0.4 for the temporal control. The parameter r in ( [To} was set to 0.8 
for all test cases. 



We run independent simulations with a variable parameter a in ( 19) for values 
a = 0.2, 0.5, 1, 2, 5, 10, 20 and 40. In Fig. [4] we plot the average long exposure 
Strehl over the field of view (red curve) and the on-axis Strehl (blue curve) for 
those a. The goal of the MCAO system is to obtain the best correction over the 
field of view, i.e., attain the largest field average Strehl. However, on-axis Strehl 
for the zenith-direction is also a quantity of interest. As can be seen from the plot, 
the peak on-axis Strehl, as well as the peak field average, is attained with a = 1 
or 2; the difference between the results can be considered negligible. 

All higher values of a over-regularize the problem and the performance of the 
algorithm, although kept stable, decreases. A low value for a corresponds to an 
under-regularized problem, and the performance of the method drops. The impor- 
tance of the regularizing term D in the presence of high photon-noise can clearly 
be observed. 



5.3. Convergence of the accelerated method. Here, we demonstrate the conver- 
gence properties of our accelerated method. We run the simulations in the same 
configuration as in the first test case, where the number of photons per subaperture 
and time-step for LGS and NGS wavefront sensors are 100 and 500, respectively. 
The readout noise is set to 3 electrons per pixel for the LGS sensors and to 5 elec- 
trons per pixel for the NGS sensors. 
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Figure 5 . Long exposure Strehl using the PCG method with dif- 
ferent number of iterations, averaged over the radial field position. 



Our accelerated algorithm is set as follows. We reconstruct three layers at the 
altitudes of the deformable mirrors using the PCG algorithm with the modified 



Jacobi preconditioner. We utilize a = 1 in ( 19 1 and choose a gain of g = 0.4 for 



the temporal control. The parameter r in ( 16 1 was set to 0.8 for all test cases. 



In Fig. [5] we plot the long exposure Strehl averaged over separation from the 
zenith after one second of simulated atmospheric propagation. We run separate 
simulations for the algorithm with 1, 2, 3, 4, 5 and 10 PCG iterations. As can 
be observed, the improvement after iteration 4 is negligible, which indicates that 
4 PCG iterations are sufficient for convergence in this configuration. 

5.4. Performance with noisy data. In this example we consider the performance 
of our methods with respect to the noise level in the LGS measurements. In other 
words, we simulate LGSs with different flux between 20 and 200. We fix the NGS 
number of photons per subaperture and time-step to 500. The readout noise is kept 
as in the previous tests at 3 and 5 electrons per pixel for the LGS and NGS sensors, 
respectively. 

We compare the performance of our methods with the matrix-vector multiply 
(MVM) method (see, e.g., [26]), which is considered to be the benchmark recon- 
structor of the ESO. The MVM is a non-iterative method in which the MAP esti- 
mate is discretized using, e.g., the Zernike polynomials. The regularized forward 
matrix is inverted and applied directly onto the measurements. The MVM that is 
presented here reconstructs three layers at DM altitudes, similar to the accelerated 
wavelet PCG method. 

We set the regularization parameter a = 1 and the gain g = 0.4. The param- 
eter r in (16 1 is tuned for each case separately O. All tests of the CG method 



are carried out with 10 iterations for the atmospheric tomography step and 4 iter- 
ations for the fitting step; the accelerated method utilizes 4 PCG iterations, which 
we found to be sufficient above. 

In Fig. [6] a comparison in quality of the reconstruction of the three methods is 
depicted. Both of the wavelet methods perform better than the MVM in almost 
all cases; the difference in the 20 photon case can be considered negligible. We 
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Figure 6. Long exposure Strehl vs. detected number of photons 
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Figure 7. Contour plots of the LE Strehl over the field of view 
(in arcsec) for MVM (left), 3-layer PCG (middle) and 9-layer CG 
(right) methods. 



believe the reason for this may be due to a better approximation of the layers using 
the wavelet basis, as well as the numerical stability of the iterative scheme, as 
opposed to matrix inversion. 

Amongst the two wavelet methods, the approach of reconstructing nine layers 
followed by a fitting step outperforms the three layer-reconstruction method in 
quality. The benefit of the full atmospheric tomography is especially emphasized 
when more photons are detected by the sensor. The disadvantage of the nine layer 
CG method is the higher computational cost over the PCG. 

To illustrate the difference between the three methods we plot Strehl values in 
the 25 directions over the field of view for MVM, 3-layer wavelet PCG and 9-layer 
wavelet CG methods for the 100 photon case in Fig. [7] 
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6. Conclusions 

In this paper we have introduced a novel reconstruction method for the atmo- 
spheric tomography problem based on wavelets. The theoretical properties of reg- 
ular wavelets enable us to apply a sparse regularization on the problem that cor- 
responds to utilizing the Kolmogorov turbulence statistics as an a priori model. 
Here, we have studied the qualitative performance of the method in the context of 
MCAO. We derived two variants of the method, concentrating more on quality by 
solving the full atmospheric reconstruction followed by a fitting step using CG or 
on speed with the layers-at-DM PCG approach. We studied the stability of the 
CG method with respect to the Bernstein-Jackson approximation. Moreover, we 
demonstrated a fast convergence of the PCG based algorithm. Lastly, we illus- 
trated the quality of the reconstructions in the low-flux regime and showed that the 
method outperforms the standard reconstruction method, called the MVM, which 
is used in the ESO simulation platform OCTOPUS. 

We believe that the wavelet method is a very promising algorithm in the field 
of atmospheric tomography. Fully utilizing the multiscale structure of wavelets 
can be approached by constructing suitable multigrid preconditioning schemes. 
Furthermore, the gain in the temporal control can be applied scale-dependently. 
Together with the careful analysis of the speed of the algorithm we leave these 
considerations for a future study. We point out that an implementation utilizing the 
discrete wavelet transform is needed in order to achieve the speed requirements of 
the E-ELT. 
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